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Abstract 

Various  theoretical  and  numerical  problems  relating  to 
helium-like  systems  in  their  ground  states  are  treated.   New 
developments  in  the  numerical  solution  of  the  Schroclinger 
equation  permit  the  solution  of  256  body  systems  with 
hard-sphere  forces.   Using  periodic  boundary  conditions, 
fluid  and  crystal  states  can  be  described;  results  for  the 
energy  and  radial  distribution  functions  are  given.   A  new 
method  of  correcting  for  low  lying  phonon  excitations  so  as 
to  extrapolate  the  energy  of  fluids  to  an  infinite  system 
is  described.   A  perturbation  theory  relating  the  properties 
of  the  system  with  pure  hard  sphere  forces  to  those  with 
smoother,  more  realistic  two  body  forces  is  introduced.  As 
in    recent  work  on  classical  systems   the  potential  is 
divided  into  two  continuous  parts;  one  is  repulsive,  one 
attractive,  the  latter  being  treated  as  a  perturbation. 
The  solution  for  the  repulsive  part  is  taken  directly  from 
the  hard  sphere  problem  when  the  radius  is  identified  as 
the  scattering  length  of  the  repulsive  part  of  the  smooth 
potential.   The  convergence  for  the  Lennard- Jones  potential 
is  very  good.   Using  our  numerical  results  for  the  hard 
sphere  problem,  with  phonon  corrections,  together  with  this 
perturbation  theory,  results  for  energy  vs.  density  agree 
with  experiment  within  our  error  of  3-10%  except  at  high 
crystal  densities.   We  carry  further  Schiff's  recent  applica- 
tion of  this  perturbation  theory  to  He   and  conclude  that 
antisymmetrization  by  the  method  of  Wu  and  Feenberg  is  the 
reason  for  lack  of  agreement  with  experiment  in  that  system. 


1.   Introduction 

In  this  paper  we  investigate  in  some  detail  the 
properties  of  the  ground  state  of  the  hard-sphere  boson 
system  and  apply  the  results  of  this  study  to  the  real 
helium  case  at  zero  temperature.   It  is  generally  assumed 
that  the  significant  features  of  the  structure  of  liquid 
helium  are  related  to  the  repulsive  part  of  the  interaction,  the 
weak  attraction  for  tne  most  part  playing  the  role  of  keeping 
the  system  bound.   A  good  qualitative  description  of  helium 
should  be  contained  in  the  hard-sphere  problem.   This  is 
our  first  motivation  to  study  that  simple  system.   In  view 
of  the  success  of  perturbation  theory  in  the  case  of 
classical  fluids  where  an  expansion  is  carried  out  around 
a  suitably  defined  hard-sphere  reference  system,  one  may 
try  to  build  the  same  kind  of  perturbation  theory  in  the 
case  of  liquid  helium,  where  the  hard-sphere  properties 
are  also  essential  inoredients.   We  shall  show  that  such 
a  theory  can  be  built  and  that  it  is  very  simple  and 
quantitatively  successful. 

The  interest  in  the  hard-sphere  Bose  system  is,  of  course, 
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not  new.   Bogolioubov,    Lee  and  Yang,    and  others, 

have  given  expressions  for  the  energy  of  the  ground  state 

which  are  exact  in  the  low  density  limit.   These  expressions 


are  unfortunately  useless  for  densities  of  physical  interest. 

The  only  existing  computations  in  the  high  density  domain 

4 
are  due  to  Hansen,  Levesque,  and  Schiff    who  use  in  the 

case  of  the  fluid  phase  a  variational  wave  function  of  the 

Bijl-Jastrow  type 


-u(r   )/2 
i     ^^  (1-1) 


:i<j 


where   u   is  a  trial  function.   The  solid  is  described  by 
multiplying  the  function  \p        by  a  product  of  one-particle 
gaussian  factors  centered  around  the  prescribed  sites  of 
the  equilibrium  lattice.   Solving  this  problem  numerically 
through  the  usual  Monte  Carlo  technique  used  in  classical 
statistical  mechanics,  HLS  obtain^   the  energy  as  a 
function  of  density  for  the  liquid  and  solid  phase,  and 
the  parameters  of  the  first-order  transition  which  lies 
between  these  phases. 

The  first  part  of  this  paper  is  devoted  to  the  solution 
of  the  quantum  hard-sphere  problem  (sections  2  through  5) . 
This  solution  based  on  a  method  devised  by  one  of  us 
is  numerical  but  essentially  exact. 

In  section  2   we  give  a  summary  of  this  numerical  method. 
It  relies  on  the  fact  that  the  wave  function  of  bosons  at 
zero  temperature,  and  the  Green's  function  of  the  Schrodinger 
equation  written  as  an  integral  equation  are  positive 
quantities  which  can  be  used  as  probabilities  in  a  Monte  Carlo 


computation.   In  the  case  of  the  hard-sphere  system  the 
only  —    but  very  great  complication   --  in  the  Green's 
function  lies  in  the  complex  boundary  conditions  imposed 
by  the   N(N-l)/2   conditions  stating  that  the  hard  cores 
do  not  overlap.   A  similar  Green's  function  can  however 
be  sampled  in  a  practical  way  on  simpler  subdomains  included 
in  the  general  complex  domain;  one  can  write  and  solve  an 
integral  equation  which  makes  it  possible  to  cover  the 
complicated  large  domain  with  the  simpler  subdomains  cmd 
thus  to  build  the  general  Green's  function.   The  use  of 
the  Jastrow  function  as  a  biasing  factor  in  the  Monte  Carlo 
method  makes  the  computation  feasible. 

Section  3   is  devoted  to  the  construction  and  sampling 
of  the  Green's  function  in  the  simple  subdomains. 

In  section  4   we  give  some  indications  of  the  technical 
aspects  of  the  method. 

The  results  of  the  computation  are  given  in  section  5. 
A  comparison   is  made  (Table  6)  with  the  variational 
results.   The  energies  obtained  through  the  "exact" 
computation  are  slightly,  although  significantly,  lower  than  the 
variational  ones.   In  the  liquid  state  the  radial  distribution 
functions    (rdf)   clearly  present  more  structure  than  the 
variational  ones.   The  off  diagonal  long  range  (ODLR) 
parameters  are  very  slightly  smaller  than  the  corresponding 
variational  ones. 

The  results  of   the  preceding  section  have  been 
calculated  for  a  small  system:   256  in  our  study.   In  section  6 


the  extension  to  quasi-infinite  systems  is  considered. 
It  requires  the  inclusion  of  long  wavelength  excitations 
which,  as  is  well  known,  play  an  essential  role  in  liquid 
helium.   Due  to  these  long  wave  length  phonons  the  energy 
of  the  dense  fluid  is  lowered  by  about  2%.   The  use  of 
the  "exact"  results  for  the  energies,  with  phonon  corrections, 
brings  a  displacement  of  the  freezing  and  melting  densities 
by  about  10%  as  compared  with  the  results  obtained  varia- 
tion ally  by  HLS. 

In  the  rest  of  the  paper,  we  apply  the  hard-sphere 
results  to  the  case  of  more  realistic  helium-like  systems. 

In  section  7  we  introduce  the  relevant  perturbation 
theory.   We  divide  the  interaction   v(r)   into  two  parts, 
as  is  done  for  the  classical  liquids:   a  part   Vq (r)  which 
gives  rise  to  the  repulsive  forces  and  which  will  be  replaced 
later  on  by  hard  spheres  located  at  the  scattering  length 
of   v„ ;   the  remainder  of  the  interaction   w(r)   which  will 
be  treated  as  a  perturbation.   The  numerical  tests  are  made 
using  the  Jastrow  wave  function  and  a  Lennard- Jones  potential 


v(r)   =   4e 
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(J)    -  (|) 


(1-2) 


for  which  many  variational  computations  have  been  made. 
The  convergence  of  the  theory  is  shown  to  be  excellent. 
Terms  of   higher  than  first  order       contribute  to  the 
energy  below  the  noise  level  of  the  computation,  i.e.  about 
0.1°  K.   The  physical  reason  for  this  good  convergence 


appears  to  be  the  following.   The  ground  state  wave  function 
of  the  system  must  be  nodeless  and  bent  as  little  as  possible 
to  keep  the  kinetic  energy  down.   This  tends  to  keep  the 
particles  as  far  away  from  one  another  as  possible.   In 
the  dense  fluid  region  this  prevents  large  density  fluctua- 
tions.  The  higher  order  terms  are  therefore  much  smaller 
than  what  they  would  be  for  a  classical  system  for  the  same 
(not  very  high)  density.   The  fact  that  the  particles  keep 
away  from  each  other  contrasts  with  the  opposite  situation 
relative  to  the  classical  system  for  which  the  radial 
distribution  function        has  a  peak  near  the  core. 
This  makes  the  quantum  case  much  simpler  than  the 
corresponding  classical  case  where  the  replacement  of  the 
repulsive  potential   Vj,{r)   by  hard  cores  cannot  be  made 
without  introducing  some  density-dependent  corrections  to 
take  into  account  the  shape  of  v^ (r) .   Here  these  corrections 
may  be  neglected;  the  density-independent  diameter  is 
obviously  the  scattering  length   a   of  the  potential   v^(r). 
To  first  order,  the  energy  of  the  total  system  is  simply 
given  as 

E  (pa^) 


of  hard  spheres  of  density   p   and  diameter   a.   The 
theoretical  justification  of  this  perturbation  theory  is 
examined  in  section  8. 


We  then  make  (section  9)   a  digression  on  the  variational 
problem  in  the  case  of  continuous  potentials:   the  perturi.a- 
tion  theory  suggests  a  form  of  the  Jastrow  wave  function 
closely  related  to  the  hard-sphere  one  and  very  different 
from  the  trial  function  usually  used  in  the  LJ  case.  We 
show  that  despite  this  difference  this  wave  function  yields 
the  usual  value,  of  the  order  of  -5.8°K   for  the  energy, 
and  an  rdf  which  presents  a  little  more  structure  than 
the  ones   i,reviously  obtained. 

In  section  10  we  give  the  results  obtained  with  the 
perturbation  theory.   We  first  show  that  using  the  varia- 
tional hard-sphere  results  obtained  with  the  HLS  trial 
function,  one  recovers  the  LJ  variational  energies  both 
for  the  liquid  and  the  not  too  dense  solid  phase.   This 
shows  the  very  large  degree  of  success  of  the  perturbation 
theory.   We  then  introduce   the  "exacf  results  for  the 
fluid   with  the  phonon  correction  included.   This  has  the 
effect  of  lowering  the  energy  of  the  liquid  by  about  0.6°  and 
of  bringing     the  location  of  the  minimum  of  the  energy 
versus  density  and  the  transition  data  into  essential 
agreement  with  experiment. 

In  the  last  section  (section  11)   we  extend  to  the 
quantum  case  the  hard-sphere  model  introduced  in  the 
classical  case  in  order  to  describe  the  structure  factor 
obtained  from  X-ray  and  neutron  scattering  experiments. 

Appendix  A  is  devoted  to  a  tabulation  of  the 
numerical  results  for  the  radial  distribution  function. 

In  Appendix  B,  we  make  some  remarks  on  the  problem 
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of  He    which  complement  the  recent  study  by  Schiff. 


2 .   Monte  Carlo  Integration  of  the  Schrodinger  Equation 


We  write  the  Schrodinger   equation  for  a  system  of 

bosons  interacting  by  hard  core  forces  in  the  following 

2 
dimensionless  form  (with   h  /2m  =  1,   potential  radius  =  a) 


-  I    ^i^(^l ^N^   =   ^*^^1 ^N^ 


(2.1) 


with  the  boundary  conditions 

\p{r,,...,r^)    =    0   when   |r^-r.|  £  a   for  all  i^j  (2.2) 

In  addition,  the  wave  function  is  periodic  in  a  box 
of  side   L  =  (N/p)"'"'^^: 


'^^(^1 ^i ^N^  =  ^(^1 VPj ^N^  ' 


p.  =  L  ><  unit  vector  in  x,y  or  z  directions 
i  =  1,.  .  .  ,N, 

Equation  (2.1)   is  written  in  the  succinct  form 


V^iij(R)   =   Eip(R)  (2.3) 


-  V^G(R,Rq)   =   6(R-Rg)  (2.4) 

Equation  (2.3)  may  be   rewritten  as 


G(R,R 


i){R)    =    E    G(R,R' )  4^(R' )  dR'  (2.5) 


If  a  succession  of  functions  is  defined  by 


^("+1)(R)   =   E  I  G(R,R')  iI^^^Nr')  dR'  ,         (2.6) 

n  =  0,1,. . . 


asymptotic  ^  for  large  n.   The  Monte  Carlo  method 

consists  in  sampling  possible  sets  of  coordinates  ^'^q^ 
at  random  from  ii         ,      an  initial  or  trial  function. 
For  any  point   R   drawn  from  ^   ^    {R   ) ,    a  point   R  ^, 
drawn  from  ip  is  obtained  by  sampling  E  G(R  ^-i  /R  ) 

considered  as  a  density  function  for   R  , ,   conditional 

-'  n+1 

upon   R  .   We  will  call  a  set  of  coordinates   R   a  "configura- 
tion" and  refer  to  the  population  of  configurations   with  a 
given  value  of  n   as  a  "generation". 

Eq  ,  the  energy  of  the  ground  state,  is  that  value  of 

(n) 
E   which  makes  the  normalization  of  i>  (i.e.,  the  size 

of  the  population  after   G   has  been  used   n   times) 

asymptotically  stable.   After  the  process  has  been  judged 

to  converge,  based  upon  estimates  of   E^   and  other  computed 

results,   E„   may  be   determined  from 

Eq   =   E-  I  4^^")  dR  /  j  4^^"^^^  dR  (2.7) 


This  process  may  be  made  computationally  more 
efficient  --  to  an  arbitrary  degree  --  by  the  following 


may  be  the  same  as  ^         (R) 


(2.8) 
Equation  (2.5)  becomes 

ip(R)  =  E'  J  [4)j{R)G{R,R')/il)j(R')]    I  (R' )    dR'  ,        (2.9) 

formally  the  same  as  (2.5).   A  sequence  of  functions 
obtained  by  successive  application  of  G  to   i|J  (R)^'    (R) 


converges  to  ^ji'ci      ^'^'^      ^n  ^^    again  that  value  that 
makes  the  population  stable.   It  will  be  shown  below  that, 

as  ip^   ^   ilj, 


I  ;p(-+iJ 


ip^"^  dR/  I  4)^"'"'  dR  (2.10) 


are  correct  on  the  average  at   n  =  0  (without  requiring 
(0) 


4)  =  tjjp.)  .  Furthermore   E.    may  in  fact  be  estimated 

with  zero  statistical  error.   One  expects  —  and  our 
experience  strongly  supports  —  that  a  reasonable 
analytical  trial  function  \p 
computational  effort  required. 

Unfortunately,  Green's  function  given  in  Eq.  (2.4) 
is  not  known  owing  to  the  complexity  of  the  boundary 
conditions.   But  all  that  is  required  is  that  a  scheme 


any   Rq .   This  is  possible  in  a  recursive  way  as  follows. 
Let   D   be  the  domain  in  configuration  space  not  excluded 


Let   Gp>(R,Rf,)   satisfy  (2.4)  but  vanish  on  the  surface 


-  V^Gq{R,Rq)  =  6(R-Rq)  (2.11) 

With  the  boundary  conditions  that  they  vanish  on  prescribed 
surfaces,  these  Green's  functions  are  symmetric.   Thus, 

-  V^G(Rj^,R)   =  6  (R-R^)  (2.12) 

Multiply  (2.11)  by  G(Rj^,R),   (2.12)  by   Gq(R,Rq)   and 
integrate  with  respect  to   R   over   D„ .   With  the  help 
of  Green's  theorem,  we  find 


The  integral  on  the  right  side  is  extended  over  the 
surface   S„  of  D„ ;   the  kernel  is  the  normal  derivative 
of   G„   which  is  everywhere  positive.  The  substitution  of 
(2.13)  into  itself  yields  a  series  of  iterated  integrals 
which  may  be  sampled  by  the  following  random  walk 
procedure . 

A  point  is  started  at   R^   and  moved  to   R'   chosen 
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then  a  new  domain   D    is  constructed,  containing   R'   anc 
a  new  point  selected  on  the  surface   S  .   The  process 
terminates  only  when  a  point  arrives  at  the  boundary  of   I 
since  (2.11)  implies 


(2.14) 


We  must  choose   D    c-.t  every  sta^e  so  that   G„  is 
known  or  can  be  sampled  and  so  that  moves  within   D 
do  not  produce  overlap  of  hard  cores.   A  particularly  con- 
venient domain  for  the  purpose  is  a  Cartesian  product 
of  N  spaces,  each  a  sphere  for  one  particle.   The  radii 
are  chosen  so  as  to  prevent  core  overlap.   This 
construction  permits   G„   to  be  sampled  in  a  straightforward 
numerical  way  (as  seen  in  section  3   below) .     As  shown 
in  reference  5,   the  requirement  that  the  wave  function 
and  hence  Green's  function  be  periodic  is  met  by  the 
computational  device  of  putting  back  on  the  left  of  the  box 
those  particles  which  leak  out  of  the  right. 

As  it  stands  (2.14)  is  satisfied,  so  in  that  form  the 
random  walk  does  not  terminate.  The  modification  required 
for  Eq.  (2.9)  removes  this  difficulty.   We  have  from  (2.13) 

^^j(\)G{R^,Rq)/^j(R^)    =    ^j(R^)Go(R^,Rq)/^j(Rq) 

+  J  {4^j(R)  [-V^Gq(R,Rq)]/iPj(Rq)  }  {i(;j(R^)G(Rj,  ,R)/^j(R)  }  dR 
^0  (2.15) 
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This  has  the  same  structure  as  (2.13)  but  a  different  kernel, 

{tJ;^[-V  G^  (R,R' )  ]/i^^}  for  which  (2.14)  does  not  hold.   In  fact, 
J    n  u         J 


-V^i|;_(R)   =   E^(R)  ijJ,(R)  (2.16) 


^0 
we  have 


•J^j^^^f-Vo^^'^o^J/'^j^^^  ^^ 


=  1-1  [Ej(R)4^j(R)  Gq(R,Rq)/iPj(Rq)  ]  dR  (2.17) 


that  the  random  walk  continues,  is  now  less  than  1. 


the  integral  on  the  right  in  (2.17)  is  that  contribution  to 


the  right  side  of  (2.9),  which  comes  from   Gq   alone. 
Thus,  when  ii        is  the  exact  eigenf  unction ,  the  process 
of  contributing  to  the  next  generation  and  that  of  going 
on  with  the  random  walk  which  develops  G  may  be  treated 
as  mutually  exclusive  random  events.   Then  the  random 
walk  always  gives  exactly  one  configuration  in  the  next 
generation  and  the  generation  size  has  no  sampling  error. 

Finally  we  observe  that  Eq.  (2.17)   is  correct  if   Gq 
is  replaced  by  G   and  the  domain  of  integration  is  D. 


Then,  since   i|^   =  0   on  the  boundary, 

r 
4|j(Rq)  =  J  Ej(R)  ^j(R)  G(R,Rg)  dR  .  (2.18) 

Again,    if      -V    ijj      =   EqI|j      , 

E^    = (2.19) 


If,  in  Eq.  (2.10)  we  set   E'  =  1   and   ij;  ^"^  =  6(R-Rq), 
then  the  expected  value  of  the  integral  in  the  denominator 
of  (2.10)  is  that  of  (2.19);  the  expected  energy  is 
identically   E„  ,   independent  of   R^ .   Taken  together 
with  the  result  stated  in  the  preceding  paragraph ,  we  see 
that  when  i)         is  the  exact  eigenf unctior\  the  energy  is 
obtained  with  neither  statistical  nor  convergence  error. 
Expectation  values  are  obtained  by  the  method  given  in 
ref.  5   in  which  the  asymptotic  size  of  the  generations 


sampling  that  led  to   R^ .   When   tjj^   is  exact,   the  estimate 
u  J 

of  the  asymptotic  size  is  carried  out  with  no  sampling  error. 
The  numerical  calculations  whose  results  are  given 


4 
by  Hansen,  Levesque  and  Schiff 


i|ij(R)  =  Tl   tanh  ((r"   -  D/b"") 
i<j 


(2.20) 


for  the  fluid,  and 
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(R)  =  exp  {-  I  A  ^  (r  -s  )^|  TT  tanh  ((r.   -  D/b""]  (2.21) 
I       i  -'  i<j  -' 


parameters   b,  m,  and  A  were  those  used  in  the  calculations 
of  ref.  4.   Computational  experience  shows  that  the  values 
which  minimize  the  energy  in  a  variational  calculation  are 
good  ones  for  our  purposes. 
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3.    Construction  and  Sampling   G^.   in  Product  Spaces 

In  the  section  above  it  is  assumed  that  at  any  stage 

in  the  random  walk,  G„(R,R  )   can  be  sampled  for  some 
0     m 

convenient  domain   D    containing   R  .   This  is  particularly 
simple  to  carry  out  if   D    is  a  Cartesian  product  of 
subspaces.   We  consider  the  case  where 

(3.1) 


and  each   d,   is  a  domain  of  the  coordinates  of  the 

k    particle  of  the  system.   The  method  holds  equally  for 

other  decompositions;  division  into  relative  and  center 

of  mass  coordinates  may  be  worth  considering  in  some  problems, 


•V^g(r  ,r'  t)  +   ^ 

and 

g(r,^,r;,t) 


Then  for   R  and  R'  in  D,  set 

Gq(R,R')   =1  TJ  g{r^,rl,t)    dt  (3.3) 


Clearly 
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and 


I'k*^}T7  9"^.'"^k'''  -" 


6  (R-R') 


(3.4) 


(3.5) 


Gq(R,R')   vanishes  for   R  on   9D   (some   r^^  on   3d^)  . 
Integrating  (3.4)  with  respect  to  time  and  using  (3.5) 
we  find  the  result  that 


{-  ^  ^k}  Gq^^'^')   =  6  (R-R-) 


(3.6) 


as  required. 

In  the  form  (3.3),  sampling   Gq   very  nearly  requires 
only  sampling   g(r,  ,rj^,t)   for  a  move  from   rj^   to   rj^ 
for  every  particle.   Consider  the  problem  of  sampling 
the  normal  derivative 


Now 


9D  =  LJ  d^  0  d2  ®  ...  ®  3  dj^  ®  ...  ®  d^^  . 


(3.7) 


Sampling  the  kernel  means  first  finding  m,  then  sampling 
points  r,  {i  ^  m)  on  the  interior  of  d^^^  and  r^  on 
the  surface  dd    .   Corresponding  to  (3.7), 


-V„G^(R,R') 


^1 


H^^m 


g(r5,r„ ,t) 


H^(t) 

~ 

H(0) 

=    1 

H(«') 

=    0 

g(r„,r„,t)  dr„ 
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(3.9) 


Green's  theorem  applied  to  (3.4)   shows  that 


dt 


3d, 


Expand       (3    .8)    into    the    form 


(3.10) 


■V^Go(R,R')   =   I[-H^(t)  T7H^(t) 
m         x-T^m 


-  h' (t) 


k^;     V^^ 


dt  . 


(3.11) 


The  sampling  procedure  is  as  follows. 

1.  Sample  t  and  m  at  random  using  the  first  quantity 
in  brackets  in  (3.11) 

2.  Conditional  on  m  and  t,  sample   r   at  random 

'^  m 

on   9d^   using  kernel   (-V^g  (r^,  r^J^,  t)  /  (-Hj^(  t)  )  . 

3.  Using  the  same  t,   sample  for  every   k  7^  m   a  point 


a  paiticle  started  at   r,   at  t  =  0 .   H,  (t)   is  the 

probability  that  it  is  absorbed  at  the  boundary  of   d, 

after  time   t   and    -  H,  (t)   is  the  rate  of  absorption  at  t. 


Suppose  for  each  k 


that  the  smallest  of  all   t,   is   t    and  that   t„   lies 


in  a  unit  interval  of  time  near   t   is 
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In  other  words  the  operations  of  selecting  the  smallest  of  k 
values   t,  ,  samples   m  and  t   as  required  for  equation  (3.11). 

Sampling   Gq(R,R')   itself  can  be  done  in  an  analogous 
way.   Of  course,  all   r,   lie  on  the  interior  of  the 
corresponding  domains   d,  .   Here  a  numerical  approximation 
to  1  r  ^k  (t)   is  used  to  sample  a  value  of  t. 

In  practice,  each  domain   d,   was  taken  to  be  a  sphere 
of  radius   aj^   centered  at   rj^ .   In  terms  of  z^   =    |r^-rj^| 
the  familiar  method  of  separation  of  variables  and  eigen- 
f unction  expansion  gives 

g(z,  ,t)  =  (2aj^Zj^)~   J  j  sin  (  tt  j  \/\)    exp  (-  tt  n  t/a^^) 
""  (3.12) 

We  see  explicitly  that  aj^  g  is  the  same  function 
of  the  reduced  coordinates  z^^/a  and  ^/\-  This  fact 
is  exploited  in  the  numerical  work. 

The  expansion  (3.12),  and  a  complementary  one  obtained  from  the 
Poisson  sum  rule  which  converges  rapidly  for  short  time, 
are  the  basis  for  various  numerical  algorithms  used  to 
sample   t,  ,  z    and  so  on.   For  the  most  part  the  z^ 
are  sampled  from  multivariate  normal  distributions  which 
can  be  done  efficiently.   Although  the  procedure  may  seem 
somewhat  elaborate,  the  computation  time  is  dominated  by 
the  usual  calculation  of  the  distances  between  particles. 
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4 .    Implementation  of  a  Computer  Program 

We  give  here  some  technical  details  of  the  computer 
program  in  which  the  ideas  outlined  in  the  preceding  section 
have  been  realized.   Readers  interested  only  in  the  theory 
may  skip  this  part  of  the  paper. 

In  any  many  body  simulation  a  crucial  detail  lies 

in  calculating  the  distances  between  pairs  of  particles. 

In  the  present  calculation  this  is  done  for  two  reasons. 

First  and  most  important,  it  is  necessary  to  know  for  each 

particle  the  distance  to  its  nearest  neighbor  so  as  to 

permit  the  choice  of  the  sphere  in  which  it  may  move  without 

core  overlap.   Second,    the  computation  of  the  value  of 

the  trial  function  \p        used  to  accelerate  the  computation 

requires  the  value  of   r. .   for  all  pairs  (cf.  eq .  (2.20)). 

For  the  latter  purpose  ii        need  not  be  as  accurate  at 

large  pair  separations  as  needed  for  good  variational 

calculations.   Hence  the  tabulation  of   r. .   is  truncated 

ID 

at  a  distance  where  the  hyperbolic  tangent  in  (3.20) 
nearly  attains  the  value  1.   The  pair  function  is  smoothly 
extrapolated  to  1  at  that  distance. 

The  tabulation  of  the  pair  distances  follows  the  proce- 
dures that  have  been  in  use  in  the  Orsay  Monte  Carlo  and 
molecular  dynamics  programs  for  some  time.   It  suffices 
to  say  here  that  tables  are  established  in  which  rather 
near  neighbors  of  each  particle  are  identified.   Most  of 
the  time  (for  a  crystal,  all  of  the  time)  distances  are 
computed  for  pairs  only  in  this  table.   When  the  tables 
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are  built  up  or  renewed  at  appropriate  intervals,  all  pairs 
are  examined. 

As  has  been  suggested  in  the  preceding  section, 
numerical  tables  are  generated  to  facilitate  sampling  at 
random  from  various  probability  distributions  required  to 
select  steps  from  the  Cartesian  product  Green's  function. 
Space  limitation  precludes  a  detailed  description  of  the 
organization  of  these  tables,  their  use,  and  the  analysis 
that  underlies  them.   We  remark  only  that  care  was  taken 
that  with  the  help  of  these  tables,  the  mappings  or  inverse 
mappings  based  on  the  Green's  function  are  very  rapid  and 
accurate.   With  respect  to  the  latter,  extensive 
numerical  tests  showed  that  in  a  step  drawn  from  the  Green's 
function  as  kernel,  the  error  in  the  mean  square  displace- 
ment is  at  most  a  few  parts  in  10  . 

The  computer   program  has  undergone  considerable  develop- 
ment during  the  period  when  it  was  being  used.   The  part 
which  has  changed  most  is  that  pertaining  to  "importance 
sampling,"  i.e.   the  transformation  of  the  dependent  variable 
from  \l)      to  \l)^      as  indicated  in  equations  (2.8)  et  seq. 
In  fact,  in  the  first  version,  this  was  not  done  although 
the  role  of  the  trial  function  in  forcing  convergence  of 
the  Green's  function  development  was  recognized  and  included 
(as  in  eq.  (2.15)).   We  estimate  that  making  the  transforma- 
tion in  a  thorough  going  way  reduced  the  computer  time 
required  by  about  a  factor  of  twenty. 


20 


In  any  case  it  becomes  necessary  to  sample  the  kernels 

(4.1) 


and 


rather  than   -V  Gj^   and   G^,   respectively.   The  techniques 
for  doing  this  have  also  evolved.   At  present  the  method 
is  based  upon  an  expansion  of   iJ^-,{R)   about   ijj  (R„).   To 
first  order 


If  we  write 


then 


tPj(R)/ti;j(RQ)   =  1  -  (R-Rq)  -Vu(Rq)  .  (4.3) 


u(R)  is  commonly  called  the  "pseudopotential" ;  we  may 

refer  to   -Vu   as  the  "pseudoforce" . 

If   ip    is  set  equal  to  a  constant  then  for  each 

sphere  the  direction  in  which  a  particle  move   is  made 

must  be  chosen  in  an  isotropic  way.   The  pseudoforce  defines  a 

favored  direction  and  a  degree  of  anisotropy  for  the 

proposed  move.   For  a  system  of  256  particles  in  a  crystal 

phase  this  elementary  use  of  the  pseudoforce  has  reduced 

the  required  computing  time  by  about  a  factor  of  10.   The 

extra  time  to  compute  the  pseudoforce  is  included.   It  is 

clear  that  further  improvements ,  possibly  resulting  in 

11 
additional  drastic  savings,  can  be  made. 
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An  awkward  feature  of  the  computational  method  lies 
in  the  necessity  to  use  a  trial  eigenvalue  and  the  role 
that  it  plays  in  the  branching  of  configurations. 
It  appears  to  require  one  to  two  hundred  generations 
for  the  energy  and  other  quantities  to  settle  down.  Since 
the  energy  changes  by  a  total  of  3  -  8%  and  since  the 
growth  of  the  population  in  n  generations  is  proportional 
to  the  n    power  of  the  trial  energy,  the  management  of 
the  calculation  in  the  face  of  statistical  uncertainties 
can  be  serious.   In  practice  runs  are  made  using  from  10 
to  80  generations;  early  in  the  relaxation  process  one 
aims  for  short  runs.   At  the  end  of  each  such  run  the 
coordinates  specifying  the  configurations  are  recorded 
on  magnetic  tape  and  used  as  input  for  the  subsequent  run. 

A  recent  series  of  runs  (crystal,  pa   =  0.4)  is 
typical.   The  process  was  followed  for  210  generations 
from  a  set  of  configurations  calculated  variationally 
after  which  an  additional  240  generations  were  followed 
during  which  the  calculated  parameters  were  constant  within 
apparent  errors  (<  1%  for  the  energy) .   A  total  of  21  hours 
of  time  on  a  CDC  6600  was  used,  about  8  in  the  early  phase. 
We  estimate  that  about  one  third  of  the  total  time  could 
have  been  saved  if  the  generation  growth  had  been  observed 
and  managed  better.   In  part,  failure  to  do  so  lies  in  the 
circumstances  under  which  the  computer  runs  were  made. 
It  is  likely  that  if  the  variance  which  is  introduced  by 
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the  branching  in  the  Green's  function  development  can  be 
further  reduced,  then  the  uncertainties  which  result 
from  the  use  of  the  trial  energy  will  also  be  reduced. 
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5.    Results  of  the  Integration  of  the  Schrodinger   Equation. 

The  numerical  results  for  the  energy  are  given  in  Table  1. 
The  errors  quoted  are  estimates  of  the  standard  derivation 
of  the  Monte  Carlo  results  of  several  runs  (3-5)  in  which 
the  properties  of  the  system  were  found  to  be  constant. 
In  all  cases  the  results  are  those  for  a  system  of  256 
particles.   For  the  fluid,  a  correction  to  an  infinite 
system  is  applied  as  explained  in  the  next  section.   Numerical 
values  of  the  radial  distribution  function  computed  for  the 
256  body  fluid  both  by  the  variational  and  the  Schrodinger 
integration  are  given  in  Appendix  A. 

The  energies  are  lov/er  tlian  those  calculated  variationally 
by  a  small  but  significant  amount.   In  the  fluid  region  the  rdf 
shows  slightly  more  structure  than  variational.   As  a  consequence 
the  main  peak  of  the  structure  factor  is  enhanced  relative  to 
the  variational  prediction. 

In  summary,  the  Jastrow  wave  function  is  found  to  be  quite 
good  as  a  first  approximation.  The  small  differences  will  turn 
out,  however,  to  be  of  decisive  importance  in  the  analysis  that 
follows . 
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6 .   Taking  the  Phonons  into  Account 

The  approach  we  have  followed  is  clearly  unable  to 
deal  with  long  wavelength  excitations ,  because  of  the 
limitation  due  to  the  size  of  the  box.   As  we  know  the 
importance  of  these  long  wavelength  excitations  in  liquid 
helium,  we  have  to  evaluate  their  effect.   In  the  framework  of  the 
Jastrow  upproximation  function   the  phonons  can  be  taken  into 
account  by  adding  to  the  short-range  pseudo-potential 
u  (r)   of  the  Jastrow  wave  function,  a  Reatto-Chester 


(6.1) 


L  2.    2^  ,  -2 

PTT  h   r  +  k 

where   c   is  the  velocity  of  sound  and   k     a  wavelength 
cutoff   which  should  be  of  the  order  of  a  few  xnterparticle 
distances   (.25/  a  is  a  guess  usually  met  in  the  literature) 
The  change  in  the  kinetic  energy  of  the  haid  spheres  is 
the  sum  of  two  terms 
2 
8m 
2 


1  =  ^ f  dr  V    u,  (r)  g^(r)  (6.2) 


,2 


^2     8i 


I  dr  (V^Ug(r)  +  V^Uj^(r))  6g(r)  .  (6.3) 


Here  the  subscript   s   refers  to  quantities  calculated 
with  the  usual  short  range  Jastrow  function.   6g(r)   is 
the  change  in  the  rdf  due  to  the  Reatto-Chester  term. 
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The  difficulty  of  this  variational  calculation  is  due 

to  the  long  range  character  of  the  Reatto-Chester  function. 

14 
It  has  been  proposed,    in  order  to  solve  this 

problem,  to  use  an  Ewald  sum  technique  as  in  the  classical 

charged  particles  system.   It  is  quite  clear  that  we  can 

thus  take  correci.ly  into  account  phonons  of  such  wave 

numbers  which  may  be  reached  with  a  finite  box  and  periodic 

boundary  conditions.   This  sets  a  limit  of  0.6  a"   in  the 

case  of  a  256  particle  system  near  solidification  density. 


effect  will  be  missed.   We  can  also  predict  the  failure 
of  the  Ewald  sum  technique  for  the  following  reason: 

in  the  fluid  region  the  variational  and  "exact"  energies 

2    2 
do  coincide  within  about  0.2  h  /ma  .   As  a  consequence 

the  contribution  of  the  short  wavelength  phonons  is  smaller 

2    2 
than   0.2  h  /ma    and  at  the  noise  level  of  the  standard 

Monte  Carlo  computation. 

In  order  to  tackle  the  problem  of  long  wavelength 

phonons,  we  shall  resort  to  a  perturbation  technique,  since 


weak  long-range  pseudo-potential  in  addition  to  the  short 


Defining 

Uj^(k)  =  J  dr  e^^'""  Uj^(r)  (6.4) 

the  dominating  terms  in  the  perturbation  expansion  are,  as 
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It  has  been  shown     that  the  proper  way  of  taking  into 
account  the  effect  of  the  size  of  the  particles  involved 
in  the  ring  was  to  put  at  each  vertex  of  the  ring  diagram 
the  structure  factor  of  the  reference  system,  with  no 
long-range  pseudopotential ,  i.e.: 


Sg(k)  =  1  +  p  j  d?  (gg{r)-l)  e^'^'''  (6.5) 


This  approach  leads  to  the  formula 


6g(r)  =  g^Cr)  (exp  (-  T  (r)  )  -  l)  (6.6) 


sf 'k)  u^  (k) 

r(k)  =  — ^ ^ .  (6.7) 

1+Uj^(k)  Sg(k) 

It  has  been  noted  by  Andersen  and  Chandler    that   r(r) 
and  thus    (5g(r)   depend  crucially  on  the  arbitrary 


repulsive  core  (where  it  has  no  physical  meaning) . 
This  dependence  appears  only  because  a  particular  subset 
of  graphs  in  the  perturbation  expansion  has   been 
selected.   AC   show  that  the  prescription   T (r)  =  0 
inside  the  repulsive  core  is  very  reasonable  and  leads 
to  excellent  results  in  the  case  of  primitive  models  of 
electrolytes  as  well  as  in  the  case  of  the  classical 
Lennard-Jones  liquids. 
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For   pa   =  0.234,   we  show  in  Figure  1  the  energy  shift 

as  a  function  of   k  .   It  is  seen  that  a  minimum  is  obtained, 
c 

for   k   =  0.3  a   ,   a  value  which  corresponds  closely  to 

the  intuitive  guess  made  above.   The  energy  is  correspond- 

2    2 
ingly  lowered  by   0.11  h  /ma  .   As  this  limiting  wave  vector 

is  substantially  smaller  than  the  smallest  wave  vector 

considered  in  the  actual  "exact"  computation  (2Tr/L  -  0.6  a 

as  we  have  seen)  we  should  also  apply  this  phonon  correction 

to  the  "exact"  energies  given  in  Table 

This  correction  on  the  energy,  although  small,  has  the 

effect  of  displacing  the  transition  parameters  of  the 

hard-sphere  system.   Let  us  consider  first  the  case  of  the 

variational  computation.    For  the  solidification  density, 

3  4     3 

one  finds   pa   =  0.24  +  0.01  instead  of    pa   =  0.23  +  0.01. 

The  melting  density  also  rises  from  0.25  +  0.01  to  0.26  +  0.01. 

We  now  apply  this  phonon  correction  calculated  in 
the  Jastrow  approximation  to  the  "exact"  results  obtained 
in  the  paper.   This  yields  the  results  obtained  in 
Table  1.    The  final  values  of  the  energy  of  the  hard 
sphere  system  are  shown  in  Figure  2  as  a  function  of  1/p. 
A  double  tangent  construction  gives  the  transition 
parameters.   We  obtain  for  the  solidification  density 
pa   =0.25+0.01,   and  for  the  melting  density 
p  a^  =  0.27  +  0.01. 
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7 .   The  Perturbation  Expansion  and  Its  Convergence 

We  shall  suppose  in  the  following  that  the  helium  atoms 
interact  through  a  two-body  potential   v(r).   This  is  not 
an  important  restriction  as  the  many-body  forces  in  helium 
are  known  to  be  weak  and  can  be  treated  as  a  perturbation. 
As  in  the  classical  theory  of  liquids    we  divide  the 
potential  into  two  parts:   a  "reference"  part   v„ (r)   and 
a  remainder   w(r)   which  will  be  considered  as  a  perturbation. 
Let   e   be  the  magnitude  of  the  potential  at  its  minimum 
and   r    be  the  corresponding  abscissa.   We  write 


=  0  for   r  >  r 


We  then  have 


w{r)  =  -  e  for   r  £  r 

=  v(r)  for   r  >  r 


(7.1) 


(7.2) 


m 


Hq  =  T  +  Vq  (7.3) 

where  Vq  is  the  potential  energy  corresponding  to  the  pair 
potential   Vq  (r)  .   Let  \pQ      be  the  normalized  wave 
function  of  the  ground  state  of  the  reference  system, 
pertaining  to  the  eigenvalue  Eq.   We  have  obviously 

<^q|Hq|4.q>  +  <4i|w|,J;>  <  <^|H|ij;>  <  <iPq|HqUq>  +  «J;q|w|^q>. 

(7.4) 
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This  can  be  rewritten  in  terms  of   g(r),  the  rdf 

for  the  total  system,  and  9n^^^       that  of  the  reference 

system, 

^  +  I  I  g(r)  w(r)  dr  <  |  1  -^  +  f  |  ^qU)    w(r)  dr  .  (7.5) 

The  r.h.s.  of  (7.5)  is  the  energy  obtained  from  first 
order  perturbation  theory.   A  complete  test  of  this 
perturbation  theory  would  require  the  knowledge  of  the 
"exact"  solution  of  the  many-body  problem  both  for  the 
reference  and  the  full  potential.   In  view  of  the  success 
of  the  Jastrow  approximation  in  the  hard-sphere  case, 
we  use  that  approximation  consistently  in  our  test  of 
the  validity  of  the  first-order  perturbation  theory. 

We  have  thus  performed,  in  the  case  of  the  Lennard- Jones 
potential  at  normal  density  of  liquid  helium,  a  series 
of  Monte  Carlo  computations  with  a  two-body  Jastrow  factor 


f(r)  =  e-(^/-)' 

We  choose  for  the  Lennard-Jones  potential  the  constants 

o  4 

e  =  10.22°,   a  =  2.556  A   appropriate  to  He  .   The  normal 

density  of  liquid  helium  corresponds  to  p  =  0.3648  a 
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Then  the  ground  state  of  the  total  system  is  obtained 
as  the  minimum  of   <T  +  V>   as  a  function  of  b,   the  ground 
state  of  the  reference  system  as  the  minimum  of   <T  +  V„>. 
In  this  last  case  the  first-order  term  is  also  obtained. 

We  thus  obtain 

E/N  =-  5.80+0.08  °K 
for 

b  =  1.17  (2)  +  0.01  a    ; 


and 


for 


and 


E/N  =  17.90  +  O.Of 


1.17  (6)  +  0.01  a    ; 


<W>/N  =  -  23.78  +  0.05    K  . 

The  errors  of  the  Monte  Carlo   computation  remain 
large  even  for  the  rather  precise  computation  reported  here. 
Within  these  errors  it  is  seen  that  the  minimum  of  <Hp.>  and 
<H>   are  obtained  for  values  of  b  which  are  very  close  as 
required  for  validity  of  the  theory. 
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Using  the  data  of  Table  2  and  the  set  of  inequalities 
(7.5)   we  can  therefore  set  an  upper  bound  of  the  order 
of  0.1°K  to  the  sum  of  higher  order  terms  in  the  perturba- 
tion theory.   We  can  also  state  that,  due  to  the  closeness 
of  the  upper  and  lower  bounds  on  the  energies,  9q^''^) 
and  g(r)   should  be  close  to  each  other. 

The  perturbation  theory  is  therefore  seen  to  work 
at  least  as  well  in  the  present  case  as  for  dense  classical 
liquids.   The  physical  reason  for  this  success  is  however 
rather  different.   In  the  classical  case  the  good  conver- 
gence of  the  perturbation  is  restricted  to  the  dense  states. 
There,  due  to  the  repulsive  cores,   the  particles  have 
little  chance  to  move  so  that  the  fluctuations  of  the 
perturbing  potential  remain  small.   In  the  quantum  fluid, 
the  density  is  relatively  small.   The  fluctuations  are 
avoided  because  of  the  tightness  of  the  wave  function: 
in  order  to  reduce  the  kinetic  energy  the  particles  must 
keep  away  from  each  other  as  much  as;  possible. 
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8.    Treatment  of  the  Reference  System 

Our  goal  is  to  replace  the  reference  system  by  a 
suitably  chosen  hard  sphere  system.   In  the  case  of  classical 
liquids  this  replacement  is  not  entirely  simple  because  the 
particles  tend  to  crowd  near  their  core,  so  that  some  account 
must  be  taken  of  the  shape  of  the  repulsion.     Because 
in  the  quantum  case  the  particles  must  stay  away  from  each 
other,  the  rdf   will  be  small  wherever  the  repulsive  potential 
is  not  zero,  so  that  all  shape  dependence  should  be 
negligible.   We  thus  anticipate  that  the  diameter   a  will 
be  determined  as  the  scattering  length   a   corresponding  to 
v„  (r)   and  that  the  energy  of  the  total  system  will  be 
given  by  the  first  order  perturbation  formula 


-^  \     dr  w(r)  g„(r/a)  ,  (8.1) 


where  E„   and  g    are  respectively  the  energy  and  the  rdf 
of  the  hard-sphere  gas  of  density  p   and  diameter  a. 

Let  us  consider  again  the  Lennard-Jones  potential 
with  parameters  given  above.   Integration  of  the 
Schrodinger  equation  yields 

a  =  0.8368  a    . 

The  normal  density  of  liquid  corresponds  to   pa   =  0.2138. 
Using  the  variational  results  of  section  5  we  obtain 
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-^  =  17.8  +  0.1°  K 


^^  =  -23.56  +  0.1°  K 

N  — 


We  therefore  obtain  for  the  total  energy  -5.76  +  0.15°  K. 
A  comparison  of  these  results  with  the  variational  results 
obtained  in  the  preceding  section  shows    excellent 
agreement. 

We  shall  now  proceed  to  give  some  theoretical  justifica- 
tion for  the  basic  expression  (8.1). 

We  shall  first  introduce  an  approximation  for  the  rdf. 

Let  us  write  in  the  hard-sphere  case 


g^(x)    =   ^l(x)    y^(x)  (8.2) 


where   i|j„  (x)   is  the  solution   of  the  two-body  hard  sphere 
problem  at  zero  energy 


r  ^  a 

(8.3) 
r  <  a 


y„(r/a)    then  is  a  smooth  function  of  r  at  the  core, 
n 

We  shall  assume  that  we  can  write 


2 

g^  (r)  =  i|;Q(r/a)  yj^(r/a)   (Assumption  A)   (8.4) 


where   i|i„  (r/a)   is  a  function  of  r  which  is  equal  to 

ijj    (r/a)       except   near   the    core.      We    anticipate    that      t|j-(r/a) 
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is  the  zero-energy  solution  of  the   two-body  Schrodinger 
equation  with  the  potential   v„(r). 

Assumption  A  is  justified  on  physical  grounds  by 
the  idea  that  only  for  close  two-body  encounters  is  g„ (r) 
different  from  the  rdf  of  a  hard-sphere  gas  with  a 
suitably  chosen  diameter. 

Assumption  A  is  not  sufficient  to  determine  the 
kinetic  energy.   In  order  to  calculate  it  we  must 
explicitly  use  the  Jastrow  approximation.   This  requires 
a  new  assumption.   We  shall  write  for  the  Jastrow  factor  fr)(r) 

fQ(r)  =  4jQ(r/a)  S^(r/a)    (Assumption  B)  (8.5) 

,     4 
where : 


S^(x)   =  f^{x)/rp^U) 


-^  tanh  - — -  (8.6) 

x-1        ,  m 


We  shall  now  see  how  these  two  assumptions  are 
connected.   It  will  turn  out  that  the  effective  expansion 
parameter  of  the  development  of  the  reference  system  around 
hard  spheres  (in  the  Jastrow  approximation  at  least)  is 

p  e^  =  p  a^  f  dx  A(x)  (8.7) 


where 


A(x)  =   ^q(x)  -  ^^(x)  (8.8) 


35 


is  negligibly  small  except  in  the  region  of  the  core. 

For  instance   at  the  normal  density  of  helium  (p„  =  0.3648  a  ) 

with  the  above  LJ  potential  it  is  found  that   pe   is 

-4 
equal  to  0.65  x  lO 

The  smallness  of  this  parameter  explains  why  we  have 

replaced  in  (8.1)  the  rdf  of  the  reference  system, 

expected  to  be  given  by  (8.4)   by  its  hard-sphere 

approximation.   For  the  same  reason.  Assumption  B 

entails  assumption  A.   Let  "Jn^^)   ^^  the  rdf  calculated 

using  assumption  B.   With  the  help  of  the  superposition 

approximation,  we  obtain 


gQ(r)    =   .J;J(|)  y^(f)  [1  +  pa^  |  dx'  g^(x')  A(x') 

X  [gjj(x-x')-l)  +  O(p^e^)  ]  .       (8.9) 
r  =  xa 
Using  the  fact  that   A(r)  is  concentrated  at  the  core, 
we  have 


pe\Q(r)  Y^U)    ^(l)  (8.10) 

where  x+1 

4)(x)   =   ^      I      dx'(g^(x')-l)  x'  .  (8.11) 


sup{ |x-i| ,o; 

t  i 
|6g^(r) 1/  gft(r)   is  smaller  than 


We  find  that,  at  the  normal  density  of  liquid  helium 
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Using  Assumptions  A  and  B,  the  energy  of  the  reference 
system  can  easily  be  cast  in  the  form 

2 

Eo  =  -^  I  dr  [\^   v'  log  f  ^  (r)  +  v^U)]    ^^(J)  ^^  (£) 


(8.12) 


de:isity   p   and  diameter   a; 
2 


-  ^^    [  dx  Vx^  log  Sjj(x)  i^li^)    -    i'li^)]  (8.13) 

4ma  ^ 


is  of  order    pe    and  can  be  safely  neglected. 


where 


(|)[ijj.]  is  seen,  a  posteriori,  to  be  sharply  concentrated 
in  the  region  of  the  core.   We  can  therefore  write 


3 

£a_ 

'2 


3        , 
|-  y^^d)    dx   $[,J;q]  .  (8.16) 


If  we  now  minimize  the  energy  with  respect  to  iJ;q  ,  we  find 
that  this  function  obeys  the  Euler  equation: 


■,,,,      -   .   .  (8.17) 
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Furthermore,  at  the  minimum,  E-   as  given  by  (8.16)  is  seen 
to  vanish.   A  direct  computation  shows  that,  at  normal 
density   E,  +  E^   is  smaller  than  .05°K. 

We  now  see  that  assumption  B  has  been  used  as  an 
intermediate  step  in  the  reasoning.   Its  explicit  use 
appears  only  in  the  evaluation  of  correction  terms  which 
are  shown  to  be  negligible.   Furthermore  it  will  be  shown 
in  the  next  section  that  assumption  B  is  very  well 
justified. 


38 


9 .   Remarks  on  the  Variational  Problem 

Assumption  B,  if  taken  literally,  has  the  advantage 
of  providing  a  Jastrow  wave  function  for  the  reference 
system  which,  given  the  variational  solution  of  the  hard-sphere 
problem,  involves  no  new  parameter.   In  order  to  put  this 
remark  in  practical  use,  we  must  give  an  analytical  form 
for  \pQ  (x)  .   Putting 

(9.1) 

we  have  found  the  following  fit  for  the  solution  of  the 
Schrodinger   equation   (8.17)   with  the  Lennard- Jones 
constants  defined  above: 

_1^ 

x^ 

(9.2) 

X  -  1  +  Cq  (1+C2(x-1)^)  exp  (-k(x-l))  for  x>l 


^0    = 

X    -    1    - 

where 

8^ 

w  = 

ioaV 

A*   = 

/mo    e 

^0    = 

0.04111 

^1   = 

0.0515 

^2    = 

-30.18 

k  =  11.56  . 
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We  can  use  (8.5)  with  (8.6)  and  (8.17)   to  perform  a 
Monte  Carlo  calculation  of  the  usual  type.   We  have  to 
compute 


(9.3) 


Using  the   Schrodinger  equation  it  is  easily  seen  that 

-2  \  i°g  ^H 


VQ(r)     h' 


2ma 
We  find,  at  normal  density 


1  (^^)^    ^  1 

2  Uq       Uq  X 


2ma 


'0 


N 


=   18.03  +  0.1 


•23.51  +  0.1 


(9.4) 


We  thus  obtain   E  =  -5.5  +  0.2  °  K.   This  should  be 
compared  with  the  value   E  =  -5.73  +  0.1   obtained 
variationally ,  and  with  the  value    E=-5.8+0.2 
obtained  from  first  order  perturbation  theory.   We 
see  that  the  agreement  is  excellent. 

In  the  past,  the  variational  computations  of  the 

energy  of  helium-like  systems  have  all  been  performed 

—  (b/r )  ^      7  8 
with  Jastrow  functions  of  the  form  e       ,   '    or 

functions  very  nearly  of  that  form.   '      All  the  functions 

hitherto  used  have  the  property  of  tending  rather  slowly 

towards  1  for  large  values  of  r.   A  Jastrow  function  of 

the  form  (8.5)  is  very  different,  and  we  wish  to  investigate 

its  use  in  calculating  the  energy  of  the  quantum  Lennard-Jones 
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fluid.   It  depends  on  one  variational  parameter  only  which 
is  the  equivalent  hard-sphere  diameter   a.   For   ij;„(r/a) 
we  shall  use  the  analytical  expression  porovided  by 
eqs .  (9.1)  and  (9.2).   With  this  Jastrow  function  we  have 
performed  Monte  Carlo  computations  for  256  particles  for 
various  values  of  the  diameter.   The  results  are  summarized 
in  Table  3. 

A  propos  of  the  energies  for  the  reference  system, 
given  in  the  second  column  of  that  table,  a  remark  is  in 
order:   In  general,  for  arbitrary  values  of  a,   the  Jastrow 
function  (9.2)  is  no  longer  a  solution  of  the  Schrodinger 
equation.   The  Laplacian  of  the  Jastrow  function  must 
therefore  be  evaluated  directly,  instead  of  using  (8.17). 
Because  (9.2)  gives  an  excellent  fit,   this  entails  no 
practical  difference. 

The  variational  computation  yields  for  the  reference 
system  an  energy  of  17.8°  K,  which  is  consistent  with  and 
not  significantly  below  the  value  17.9°  K  obtained  with 
the  function  exp[-(b/r)  ].   The  energy  minimum  is  obtained 
for   a  =  0.82  (5).   This  is  close  to  (but  a  little  smaller 
than)  the  scattering  length,  especially  if  one  realizes 
that  at  the  minimum  the  energy  varies  but  little  with  a. 
Assumption  B  of  the  preceding  section  is  well  verified. 

The  total  energy  per  particle  is  equal  to  -5.82  +  0.1  "  K. 
This  turns  out  to  be  little  different  from  the  energy  obtained 
with  the  Jastrow  wave  function  exp[-(b/r)  ],  i.e.  -5.73^0.1  °    K, 
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Taking  into  account  the  usual  tendency  to  underestimate  the 

errors  in  Monte  Carlo  computations,  we  do  not  find  convincing 

evidence  in  the  literature  that  a  lower  value  of  the  energy 

can  be  obtained  with  a  Jastrow  wave  function  (at  least  when 

long  wave  length  phonons  are  not  included).   The  form  (9.2) 

_  ,,  ,  ,  m 
for  the  Jastrow  function  is  sounder  than  the  usual  e    ^''^ 

because  it  exhibits  hard  sphere  features  which  are  the 

most  important  characteristics  of  liquid  helium.   We  compare 

in  Figure  3  the  structure  factors  obtained  with  these 

two  different  Jastrow  factors.   The  overall  agreement 

is  good.   As  expected,  with  a  hard-sphere-like  Jastrow 

factor,  one  gets  a  higher  main  peak  in  the  structure 

factor  than  with  a  smoother  trial  function. 
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10.   Results   of  the  Perturbation  Theory 

We  shall  apply  perturbation  theory  to  the  case  of 
the  Lennard- Jones  potential.   As  in  that  case  complete 
and  reliable  results  exist  only  for  short  range  Jastrow 
functions,  we  shall  as  a  first  step  use,  in  order  to  be 
consistent,  the  corresponding  variational  hard-sphere 
results  neglecting  phonons  in  both  cases.   The  corrections 
will  be  made  afterwards. 

We  use  the  values  of  the  variational  rdf  listed 
in  Appendix  A    in  order  to  calculate,  with  the  help  of 
formula  (8.1),  the  energy  of  the  quantum  fluid  composed 
of  Lennard-Jones  molecules.   This  is  compared  in  Table  4 
with  the  variational  results  of  ref.  8   for  the  liquid 
and  with  those  of  Hansen  and  Pollock    for  the  solid. 
The  statistical  errors  on  the  energy  and  rdf  of  the  hard-sphere 
gas  entails  an  error  on  the  energy  obtained  from  perturbation 
theory,  which  is  difficult  to  estimate  because  there  is 
clearly  a  correlation  between  the  errors  on  the  zeroth  and 
first  order  terms.   For   pa  =  0.2138  we  have  made  two 
independent  runs ,  the  results  of  which  are  compared  in 
Table  4.   The  final  energies  differ  by  0.1  °    K.   The  error 
in  the  solid  region  is  larger  because  there  we  replaced 
the  rdf  for  distances  larger  than  2.4  a   by  1 ;  the  effect 
of  this  cut-off  is  clearly  not  negligible. 

We  expect  the  perturbation  theory  to  break  down  for 
low  densities  where  one  may  get  large  density  fluctuations 
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without  increasing  the  kinetic  energy  much  and  in  the  high 
density  solid  region  where  the  hard-sphere  approximation 

must  break  down.   There  is,  in  contrast  with  the  classical 

21 
case,    a  large  domain  of  density  where  the  perturbation 

theory  with  the  hard-sphere  approximation  works  well  in  the 

case  of  the  solid.   The  reason  for  this  success  is  that 

in  the  quantum  case  solidification  occurs  at  a  quite  low 

density  because  the  particles  must  stay  away  from  the  core 

of  their  neighbors,  and  the  details  of  the  repulsion  are 

not  seen.   As  a  whole,  we  do  not  see  in  Table  4   any  clear 

discrepancy  between  the  perturbation  results  and  the  results 

obtained  variationally .   There  may  be  an  exception  to  that 

statement  for  the  lowest  liquid  density   (p  =  0.167)  where 

the  rather  large  discrepancy  may  be  due  to  our  neglect  of 

density  fluctuations,  and  for  the  highest  solid  density  where 

the  large  discrepancy  between  the  perturbation  and  variational 

results  may  be  significant.    As  a  whole  it  appears  that 

the  perturbation  theory  applies  in  a  density  domain  which 

is  remarkably  large. 

Now  that  the  validity  of  the  perturbation  theory  has 

been  checked  within  the  framework  of  the  Jastrow  approximation, 

we  shall  replace  the  approximate  HS  results  by  "exact"  ones. 

There  are  three  kinds  of  corrections  which  all  tend  to  lower 

the  energy. 
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1)  The  replacement  of  the  variational  HS  energies 
by  "exact"  ones,  as  computed  for  256  particle  systems. 
This  appears  to  be  the  leading  correction. 

2)  The  use  of  "exact"       instead  of  variational  rdf 
in  the  computation   of  the  first-order  term.   This  correction 
is  negligible  in  the  solid  region  and  very  small  (of  the 
order  of  0.05°)  in  the  fluid  region. 

3)  The  phonon  correction  is  small  but  not  completely 
negligible.   With  the  value  of  the  HS  diameter  which  has 
been  chosen,  and  with  the  Michels   constants  for  the  LJ 
potential,  we  obtain  for  the  velocity  of  sound  in  the 
hard-sphere  system   the  value  295  m/s  for  the  equilibrium 
density  of  liquid  helium.   This  is  not  far  from  the 
experimental  value   (2  39  m/s) .   Due  to  the  small  size  of 
the  phonon  correction,  it  is  rather  irrelevant  which 
value  we  use. 

The  perturbation  results  thus  obtained  are  shown  in 

22  23 
Table  5  and  compared  with  the  experimental  values.   ' 

It  is  seen  that,  within  the  estimated  error,  which  is  by 

no  means  small,  there  is  no  clear  discrepancy  except  at 

the  highest  solid  density  where  the  Lennard-Jones  potential 

clearly  gives  energies  which  are  too  low  as  already  known 

from  the  work  of  Hansen  and  Pollock.     We  can  obtain  a 

confirmation  of  these  results  in  the  following  way:  as  v/e 

have  seen,  the  correction  to  be  applied  to  the  variational 

results  is  essentially  that  in  the  kinetic  energy.  \}e    can 
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therefore  estimate  the  exact  energy  by  adding  to  the  variational 
LJ  energies  a  term 

6E  =  exact  energy  of  the  equivalent  HS  system  minus 
variational  energy  of  the  HS  system. 

In  this  process,  we  use  smoother  data  than  in  the 
straight  perturbation  theory.   The  LJ  variational  data 
are  quite  precise  and  smooth,  and   6E  can  also  be  smoothed. 
The  values  of   6e  which  are  used  are  given  in  Table  6 
and  the  energies  thus  obtained  are  given  in  column  5 
of  Table  5.   They  are  entirely  compatible  with  the  perturba- 
tion energies  and  with  experiment  (except  at  very  high 
density) . 

We  did  not  try  to  determine  from  the  above  data  the 
transition  and  equilibrium  densities.   They  are  obviously 
compatible  with  the  experimental  data,  within  a  large 
uncertainty  which  may  amount  to  10%. 
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11.   The  Hard  Sphere  Model 

In  the  preceding  sections ,  the  close  relationship 
between  the  reference  system  and  the  hard-sphere  fluid 
has  been  exhibited.   It  has  moreover  been  shown  that  the 
attractive  forces  change  the  structure  only  a  little. 
In  particular,  it  has  been  demonstrated  in  section  9, 
on  the  basis  of  variational  computations,  that  a  hard-sphere- 
like  Jastrow  function  is  perfectly  acceptable  for  the 
description  of  the  full  L.T  fluid.   Still  within  the 
variational  framework,  we  can  obtain  a  further  test  of 
the  validity  of  this  hard  sphere  representation  by 
comparing  the  off-diagonal  long-range  order  parameter 
obtained  '   for  the  LJ  fluid  at  normal  density  (with  the 
^-(b/r)   Jastrow  factor),  i.e.   n^  =  0.105  +  0.005, 
with  the  HS  variational  value  at   pa   =  0.213  8  where  we 

These  values  are  clearly 

compatible.   Last,    we  may  try  the  solidification 

24 
criterion  used  in  the  classical  case:     the  real  system 

solidifies  when  the  equivalent  HS  system  does.   The 

solidification  density  of  the  HS  system  in  the  variational 

approximation  is  obtained  for   pa   =0.23+0.01.   This 

gives  for  the  LJ  system   pa   =0.39+0.2,   which  should 

be  compared  with  the  value   pa   =  0,38   obtained  through 

the  variational  computation  made    directly 

in  the  LJ  system. 
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Let  us  now  go  beyond  the  Jastrow  approximation. 
In  Figure  4  we  show  the  "exact"   S (k)   appropriate  to 
the  normal  density  of  liquid  helium.   We  take  for  the 
full  LJ  system  the  same  equivalent  hard-sphere  density 
pa   =  0,2138   corresponding  to  a  diameter  a  =  0.8368  a 
as  in  the  case  of  the  reference  system.   This  choice  is 
motivated  as  follows:   we  know  from  the  computations  of 
section  9   that  the  attractive  forces  have  the  tendency 
to  reduce  the  equivalent  diameter  to  a  value  around 
a  =  0.81  a.   On  the  other  hand,  we  know  that  in  order 
to  describe  real  helium  more  realistic  potentials  than 

the  usual  LJ  potential  used  in  this  study  should  be  used 

25 
such  as  the  potential  due  to  Beck     or  the  modified 

Lennard- Jones  potential  (LJ2)  considered  by  Hansen 

and  Pollock.     These  potentials  have  a  wider  core  than 

the  LJ  potential  with  Michels  constants.  They  both  yield  the 

value  a  =  O.860      as  the  equivalent  hard-sphere  diameter 

for  the  reference  fluid.   We  thus  see  that  keeping  the 

value   a  =  0.8368a    realized  a  compromise  between  the 

above  two  effects. 

The  determination  of  the  "exact"   S (k)   involves  some 

manipulations  which  should  be  indicated.   For   pa   =0.2, 

we  can  compute  precisely  the  difference  between  the  exact  and 

variational  rdf.    This  difference  is  scaled  and  added  to 

the  variational  rdf  at   pa   =  0.2138.   In  order  to  get 

the  Fourier  transform  of  the  rdf,  we  have  obtained  its  value 

beyond  r  =  2.4a    by  using  the  Ornstein-Zernike  relation 
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and  the  assumption  that,  for  large   r,   the  direct  correla- 
tion function  is  negligible.    The  phonon  contribution  is 
then  calculated  as  indicated  in  section  6,  using  for 
the  sound  velocity  the  experimental  value  and  for  the  cut-off 

wave  vector  the  value   k   =  0.3  a     that  leads  to  the 
c 

energy  minimum. 

The  structure  factor  thus  determined  is  compared 
in  Figure  4  to  experimental  results  obtained  by  X-ray 

scattering.   Hallock's  data    have  been  used  for 

o_-|  27 

k  <  1.1  A      (ka  <  2.3).   Achter  and  Meyer's  results 

are  shown  for  higher  values  of   k.   It  is  seen  that,  at 

low   k    our  results  agree  well  with  those  of  Hallock 

although  they  show  the  shoulder  predicted  by  Miller,  Pines 

2  8 
and  Nozieres     a  little  more.    We  may  note  that  this 

shoulder  was  not  observed  by  Achter  and  Meyer  whose  results 

may  differ  from  Hallock's  by  as  much  as  10%,   In  the  k 

region  corresponding  to  the  rise  of  the  first  peak,  the 

experimental  and  theoretical  results  agree  quite  well. 

The  position  of  the  first  peak  and  the  zeros   in  S(k)-1 

differ  markedly.   We  note  however  that  this  disagreement 

is  no  longer  present  when  one  considers  the  results 

29 
obtained  from  neutron  scattering  by  Henshaw.     There  however 

the  first  peak  rises  to  a  much  higher  value  than  in  Figure  4. 

In  order  to  get  more  insight  into  the  origin  of  the 

discrepancy  at  large  k,  we  reexamine   the  experimental  data 

used  in  Figure  4.   The  rdf  obtained  from  these  results  by 
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Fourier  transform  is  not  strictly  zero  inside  the  core. 
If  we  correct  this  rdf   by  imposing  simply  that  it   vanish 
for   r  ^  a,   the  structure  factor  that  we  obtain  by  Fourier 
transform  is  no  longer  identical  with  the  experimental  one. 
We  can  impose  a  compatibility  with  experiment  in  the  low 
k  region  (including  the  first  peak  and  a  strict  vanishing 
of  the  rdf  inside  the  core  if  we  iterate  the  procedure 
described  above,  based  on  a  succession  of  Fourier  transforms. 
In  this  way,  a  self-consistent  structure  factor  is  recon- 
structed, with  admittedly  some  measure  of  arbitrariness. 
The  results  are  shown  in  Figure  4.   It  is  seen  that  the 
agreement  with  the  hard  sphere  structure  factor  is  somewhat 
improved.   We  believe  the  remaining  discrepancy  to  be  due 
more  to  the  errors  in  the   hard-sphere  rdf  and  the  experi- 
mental structure  factor  than  to  a  defect  in  the  model. 

The  "exact"  ODLRO  parameter  is  estimated  as 
nQ  =  0.095  +  .001  at  pa   =0.2.   It  is  very  different  from 
the  recent  estimate    of  n   =  0.0241  +  .01   which  is 
obtained  indirectly  from  experiment  at  1.2  °  K. 

Last,  we  examine  the  solidification  criterion.   The 
freezing  density  for  hard  spheres   pa   =0.25+0.01   yields 
the  value   pa   =  0.427.   At  the  solidification  density, 
the  main  peak  of  the  hard  sphere  structure  factor  reaches 

the  value  1.40  which  should  be  compared  to  the  value  2.85 

24 
in  the  classical  case.     A  value  of  the  order  of  1.40  should 

therefore  be  the  maximum  value  reached  by  the  structure 

factor  of  liquid  hf.>lium  at   very  low  temperature  in  the  whole 

fluid  range. 
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Appendix  A 

We  tabulate  the  radial  distribution  function  for 
the  hard  sphere  Bose  fluid  as  computed  by  the  Monte  Carlo 

integration  of  the  Schrodinger   equation  and  by  the  vari- 

4 
ational    method  with  the  Jastrow  factor 


Results  are  given  at   pa   =  0.2,  0.244,  and  0.27. 
The  numbers  tabulated  as   g   are  in  fact  average 
values  over  the  interval  ending  in  the  r.   That  is 


n  n 

g(r^)  =   I    r2g(r)  dr  /   | 


We  estimate  the  statistical  error  of  the  values  as 
less  than  2%   except  for  the  first  few  entries  which 
have  errors  of  5-10%.   A  systematic  error  of  a  few 
percent   at  the  peak  of   g(r)   is  possible  owing  ^:o 
incomplete  convergence  of  the  iteration  process 
used  in  the  calculation. 
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Appendix  B 

The  problem  of  solid   He   is  treated  exactly  as  above. 

As  long  as  we   ne/lect  the  (small)  exchange  effects,  the 

4  3 

only  difference  with   He   lies  in  the  lighter  mass  of  He  . 


reference  part  of  the  LJ   potential  in  the  He   case  and  found 
the  value  0.82  3  a.   He  has  shown  that,  in  the  variational 
framework,  the  perturbation  theory  works  as  well  as  in  the 


with  "exact"  results  for  the  HS  boson  solid  instead  of 
variational   results.    As  above,  we  give  both  the 
results   obtained  by  a  straight  application  of 
perturbation  theory  and  the  estimate  obtained  from 
correcting  the  variational  results  for  the  difference 
between  "exact"  and  variational  HS  results.   Both  methods 

lead  to  energies  which  agree,  within  the  rather  large 

4 
errors,  with  experiment,  as  in  the  He   case. 

For  the  liquid  state,  the  antisymmetry  of  the  wave 

function  must  be  taken  into  account.   In  the  absence  of 

"exact"  hard  sphere  results  in  the  fermion  case,  one  is 

led   to  rely  on  the  Wu-Feenberg    approximation: 


He^"^^ 


where  ^„      is  the  solution  of  the  mass  3  boson  problem 
and  i>        is  the  appropriate  Slater  determinant  of  plane 
waves.   As  Wu  and  Feenberg   have  shown,  if  the  ground  state 

52 


properties  of  the  boson  system  are  known,  the  energy  of 

the  He   fluid  can  be  calculated  by  making  a  cluster  expansion 

of  the  determinant.   This  expansion  seems  to  converge  quite 

well.   For  the  Lennard- Jones  potential  variational  results 

9 

are  available.     In  that  case,  Schiff  has  shown 

that  the       perturbation  methods  also  work  quite  well. 
Schiff 's  results  rely  on  three  assumptions:  (1)  variational 
results  are  used  for  the  hard-sphere  boson  fluid. 

(2)  The  Lennard- Jones   potential  is  a  good  approximation 
for  the  He  -  He  interaction. 

(3)  The  Wu-Feenberg  ansatz   is  a  good  approximation. 

Here  (Table  8)   we  correct  for  (1)  by  using  "exact" 

4 
hard  sphere  results.   We  know  from  the  He   study  that 

the  error  entailed  by  the  use  of  the  LJ   potential  is 

small  (less  than  0.5°).   The  very  clear  discrepancy 

revealed  in  Table  8    in  the  He    is  then  largely  due 

to  the  use  of  the  Wu-Feenberg  ansatz,  whose  effect 

amounts  to  an  overestimation  of  the  energy  by 

at  least  1  ° . 
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state 

pa^ 

^1 

phonon 
correction 

^0 

fluid 

0.166 

4.24  + 

.05 

-  0.09 

4.15  + 

.05 

0.200 

5.80  + 

.05 

-  0.11 

5.69  + 

.05 

0.244 

8.28  + 

.09 

-  0.14 

8.14  + 

.09 

0.270 

10.65  + 

.07 

-  0.15 

10.50  + 

.07 

F.C.C. 

0.244 

8.50  + 

.05 

crystal 

0.270 
0.300 
0.400 
0.500 

10.12  + 
12.27  + 
21.26  + 
34.45  + 

.04 

.07 

.07 

.25 

1 

2         2 

Energies    in    units    of      ft    /ma         for   the   ground   state 

of  the  quantum  hand  sphere  fluid  and  F.C.C.  solid. 


integration  of  the  Schrodinger  equation  for  a  system 
of  256  particles.   The  phonon  correction  to  an  infinite 
system  is  obtained  by  the  method  of  section  6. 
For  the  crystal  the  results  for  256  particles  are  given 
with  no  corrections.   Errors  are  estimates  of  the 
standard  deviation. 
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Table  2 


b 

(<T>/N)  ° 

f:VQ>/N)° 

(Eq/N)° 

(<W>/N)  ° 

(E/N)  ° 

n 

1.16 

13.59 

4.45 

18.04 

-23.83 

-  5.79 

650000 

1.17 

13.95 

4.09 

18.04 

-23.80 

-  5.76 

650000 

1.18 

14.30 

3.60 

17.90 

-23.76 

-  5.86 

400000 

1.20 

15.07 

3.06 

18.15 

-23.68 

-  5.53 

650000 

Energies  obtained  variationally   for  the 
reference  system  (column  4)   and  for  the  system 
with  the  full  LJ   potential  (column  6) . 
n   is  the  total  number  of  configurations. 
The  statistical  error  on  the  kinetic  and 
potential  energies  is  of  the  order  of  0.05  °. 
The  error  on  the  total  energies  is     about  O.OE 
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Table  3 


a 

(Eo/N)° 

(  E/N  )  ° 

n 

0.75 

20.64  +  0.08 

-  3.64  +  0.12 

600  000 

0.80 

18.04  +  0.04 

-  5. 78  +  0.06 

2  500  000 

0.82 

17.84  +  0.04 

-  5.81  +  0.06 

2  300  000 

0.8368 

18.04  +  0.04 

-  5.48  +  0.06 

2  000  000 

0.86 

18.40  +  0.04 

-  4.93  +  0.06 

2  300  000 

Variational  energies  for  the  LJ  system 

using  the  Jastrow   factor  defined  by  (9.1)  and  (9.2) 

Column  1:   equivalent  hard  sphere  diameter. 

Column  2  and  3:   energy  of  the  reference  system 

and  of  the  complete  system   respectively. 
Column  4:   total  number  of  configurations  in  the 

Monte  Carlo  computation. 
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Table  4 


1 

3 

pa 

(E„/^° 

(cW>/M  ° 

(E/N  )  ° 

^ar/^^° 

0.167 

0.283 

4.32+0.1 

11.38 

-16.64 

-    5.26+0.2 

-5.67+0.1 

0.2 

0.341 

6.01+0.1 

15.83 

-22.09 

-    6.16+0.2 

-5.91+0.1 

0.2138 

0.364 

6.85+0.1 

18.05 

-23.89 

-    5.84+0.2 

-5.73+0.1 

0.2138 

0.364 

6.75+0.1 

17.80 

-23.56 

-    5.76+0.2 

-5.73+0.1 

0.234 

0.4 

7.92+0.1 

20.87 

-26.27 

-    5.40+0.2 

-5.25+0.1 

0.244 

0.416 

8.86+0.1^23.33 

-28.44 

-    5.11+0.4 

-4.97+0.2 

0.27 

0.461 

10.6   +0.1 

27.92 

-32.73 

-    4.79+0.5 

-4.75+0.2 

0.3 

0.512 

12.69+0.1 

33.40 

-36.50 

-    3.10+0.6 

-3.87+0.3 

0.35 

0.597 

16.8   +0.2 

44.20 

-44.56 

-    0.36+0.8 

-1.4    +0.4 

0.4 

0.683 

21.9    +0.2 

57.70 

-54.85 

2.87+1 

3.2    +0.5 

0.5 

0.853 

35.5   +0.3 

93.5 

-66.02 

27.5    +2 

22.8    +    1 

Energies  of  the  LJ   system  computed  within 
the  Jastrow  approximation.   Column  6  gives 
the  result  of  the  perturbation  theory. 
Column  7  gives  the  variational  results. 
The  upper  half  of  the  table  refers  to  fluid 
states,  the  lower  half  to  solid  states. 
For  the  density   pa"^  =  0.2138,  the  results 
of  two  independent  Monte  Carlo  runs  have  been 
given  for  comparison. 
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Table  5 


pa^ 

Pa^ 

(^var/N)° 

(Vrt/^^° 

(E^^^/N)° 

■ 

0.167 

0.283 

-    5.67+0.1 

-    5.72+0.2 

-    6.13+0.2 

0.2 

0.341 

-    5.91+0.1 

-    7.14+0.2 

-    6.77+0.2 

(-    7.09) 

0.244 

0.416 

-    5.02+0.1 

-    6.50+0.2 

-    6.67+0.2 

-    6.84 

0.244 

0.416 

-    4.97+0.2 

-    6.11+0.4 

-    5.97+0.4 

(-    5.5    ) 

0.27 

0.461 

-    4.75+0.2 

-    5.79+0.5 

-    5.75+0.5 

-    6.08 

0.3 

0.512 

-    3.87+0.3 

-    4.10+0.6 

-    5.00+0.6 

- 

-    5.3 

0.4 

0.683 

3.2    +0.5 

1.2+1 

1.5    +1 

2.1 

0.5 

0.853 

22.8    +1 

24.8    +2 

20.1    +2 

26 

Results  of  perturbation  theory.   Column  4: 
perturbation  energies  using  hard  sphere  results 
Column  5:   Estimated  energies  obtained  by  adding 
to  the  variational  LJ   results  of  column  3  the 
hard  sphere  correction   6E   given  in  Table  6. 
The  upper  half  of  the  table  refers  to  fluid  states 
the  lower  half  to  solid  states. 
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Table  6 


pa^ 

^exact/N 

^a/^ 

6E/N 

0.166 

0.2 

0.244 

4.15 
5.69 
8.1J4 

4.32 
6.01 
8.80 

0.17 
0.32 
0.66 

0.244 

0.27 

0.3 

0.4 

0.5 

8.50 
10.12 
12.27 
21.26 
34.45 

8.86 
10.50 
12.65 
21.9 
35.5 

0.38 
0.38 
0.43 
0.64 
1.05 

Smoothed  difference   5E/N   between  the  "exact" 
energy  per  particle  of  the  hard  sphere  system 
and  its  variational  estimate.  Energies  in  units 
of   h  /ma  . 
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1      ^ 

5s 

^v 

r 

^s 

^v           ^ 

^s 

^v 

r 

% 

^v 

(1.10 

0.049 

0.044 

2.20 

1.025 

1.012    3.30 

1.045 

1.021 

4.40 

1.016 

0.995 

11.20 

0.239 

0.259 

2.30 

0.949 

0.969    3.40 

1.037 

1.018 

4.50 

0.993 

0.997 

1.30 

0.542 

0.567 

2.40 

0.912 

0.94li  3.50 

1.009 

1.017 

4.60 

1.000 

1.002 

1.40 

0.870 

0.876 

2.50 

0.880 

0.9331 3.60 

0.990 

1.015 

4.70 

0.990 

0.998 

1.50 

1.113 

1.121 

2.60 

0.  890 

0.947! 3.70 

0.990 

1.005 

4.80 

0.999 

1.004 

1.60 

1.280 

1.253 

2.70 

0.907 

0.943| 3.80 

1.005 

0.998 

4.90 

0.999 

1.001 

1.70 

1.394 

1.299 

2.80 

0.958 

0.964    3.90 

0.992 

0.997 

5.00 

1.007 

1.001 

1.80 

1.322 

1.274 

2.90 

0.990 

0.981    4.00 

0.993 

0.999 

5.10 

0.987 

1.000 

1.90 

1.236 

1.198 

3.00 

1.025 

0.992    4.10 

0.995 

0.993 

5.20 

1.003 

1.002 

2.00 

1.198 

1.142 

3.10 

1.045 

1.005    4.20 

0.990 

0.989 

5.30 

1.005 

1.002 

2.10 

1.065 

1.069    3.20'l.04  8ll.01l"4.30'0.998'0.996'5.40 

1.00311.000 

pa^    =    0.244 

1.10 

0.044 

0.031 

2.20 

0.912 

0.923 

3.30 

1.119 

1.103 

4.40 

0.990 

0.988: 

1.20 

0.252 

0.245 

2.30 

0.820 

0.847 

3.40 

1.077 

1.073 

4.50 

1.014 

1.012 

1.30 

0.587 

0.636 

2.40 

0.78  6 

0.801 

3.50 

1.044 

1.042 

4.60 

1.048 

1.040 

1.40 

0.955 

1.035 

2.50 

0.787 

0.790 

3.60 

0.998 

1.004 

4.  70 

1.057 

1.059 

1.50 

1.271 

1.285 

2.60 

0.812 

0.820 

3.70 

0.963 

0.976 

4.80 

1.066 

1.060 

1.60 

1.457 

1.414 

2.70 

0.880 

0.881 

3.80 

0.936 

0.952 

4.90 

1.056 

1.057' 

1.70 

1.504 

1.428 

2.80 

0.952 

0.953 

3.90 

0.931 

0.940 

5.00 

1.036 

1.035! 

1.80 

1.409 

1.365 

2.90 

1.021 

1.018 

4.00 

0.933 

0.935 

5.08 

1.021 

1.023 

1.90 

1.283 

1.291 

3.00 

1.085 

1.077 

4.10 

0.943 

0.937 

2.00 

1.169 

1.162 

3.10 

1.129 

1.110 

4.20 

0.949 

0.952 

2.10 

1.029 

1.037 

3.20 

1.131 

1.125 

4.30 

0.969 

0.966 

pa"*    =    0.27 

1.05 

0.017 

0.021 

2.05 

0.919 

0.938 

3.05|1.054 

1.035 

4.05 

1.005    1.003 

1.10 

0.100 

0.109 

2.10 

0.  886 

0.909 

3.10 

1.055 

1.026 

4.10 

1.008    1.000 

1.15 

0.256 

0.280 

2.15 

0.877 

0.889 

3.15 

1.036 

1.024 

4.15 

1.024    1.008 

1.20 

0.450 

0.507 

2.20 

0.874 

0.877 

3.20 

1.  004 

1.014 

4.20 

1.0061. 004 

1.25 

0.686 

0.739 

2.25 

0.843 

0.878 

3.25 

!l.024 

1.004 

4.25 

1.014 ;0.998 

1.30 

0.931 

0.990 

2.30 

0.863 

0.  887 

3.30 

.0.994 

1.001 

4.30 

1.017:1.005 

1.35 

1.154 

1.184 

2.35 

0.877 

0.894 

3.35 

:0.994 

1.011 

4.35 

1.00311.005 

1.40 

1.  320 

1.330 

2.40 

0.887 

0.910 

3.40 

10.999 

0.997 

4.40 

0.996!l.002 

1.45 

1.415 

1.425 

2.45 

0.924 

0.928 

3.45    1,017 

0.981 

4.45 

O.993I1.OO9 

1.50 

1.489 

1.469 

2.50 

0.922 

0.941 

3.50    0.997 

0.983 

4.50 

1.00811.003 

1.55 

1,527 

1.492 

2.55 

0.942 

0    968 

3.55    0.986 

0.985 

4.55 

0.996ll.003 

1.60 

1.478 

1.441 

2.60 

0.961 

0.972 

3.60 

0.972 

0.936 

4.60 

1.000|0.999 

1.65 

1.431 

1.387 

2.65 

0.987 

1.000 

3.65 

0.963 

0.9  84 

4.65 

0.999'0.994 

1.70 

1.364 

1.313 

2.70 

1.022 

1.005 

3.70 

0,975 

0.990 

4.70 

I.OIC    1.002 

1.75 

1.273 

1.243 

2.75 

1.034 

1.022 

3.75 

0.970 

0.985 

4.75 

0.995    0.998 

1.80 

1.217 

1.162 

2.80 

1.060 

1.022 

3.80 

0.976 

0  .  9  89 

4.  80 

1.006(0.999 

1.85 

1.108 

1.094 

2.85 

1.034 

1.044 

3.85 

0.988 

1.001 

4.85 

I.OO2I0.999 

1.90 

1.049 

1.025 

2.90 

1.062 

1.050 

3.90 

0.984 

0.993 

4.90 

0.99710.996 

1.95 

0.986 

1.004 

2.95 

1.033il.053 

3.95 

0.992 

0.997 

4.91 

1.008    0.995 

2.00 

0.957 

0.965 

3.00jl.033jl.035 

4.00 

1.001 

0.994 

Radial  distribution  function  for  fluid  states  of  the  hard  sphere 
boson  system  calculated  by  variational  method  (g^)  and 
Schrodinger   integration  (g  ) • 
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3 

pa 

PO^ 

^^var/^° 

'■=pert/"'° 

(E^^^/N)o 

E         /N)  ° 

0.167 

0.283 

-1 

-1.2 

-1.3 

-    2.4 

0.2 

0.341 

0 

-0.35 

-0.75 

-    2.2 

0.2138 

0.364 

0.6 

-0.4 

-    1.8 

0.234 

0.4 

2.2 

0.75 

0.244 

0.416 

3.2 

0.55 

1.65 

0.244 

0.416 

1.5 

0.2   +0.6 

0.1    +0.6 

-    0.4 

0.27 

0.461 

3.2 

1.8   +0.8 

1.8   +0.8 

1.1 

0.3 

0.512 

5.8 

4.2    +1.0 

4.3   +1.0 

4.6 

0.35 

0.597 

11.7 

9.8   +1.2 

11.5 

Results  of  the  perturbation  theory  applied  to  He  . 
Same  as  in  Table  5.   The  results  for  the  liquid  state 
have  been  obtained  using  the  Wu-Feenberg  ansatz. 
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Figure  Captions 


Fig.  1     Energy  shift  due  to  long  wavelength  phonons  in 

2    2 

units   h  /ma   ,   as  a  function  of  the  cut-off 

wave  vector   k    for  the  density   pa   =  0.234. 


Fig.  2    "Exact"  energy  per  particle  as  a  function  of  1/pa  . 


Fig.  3     Structure  factor  for  the  Lennard- Jones  system  as 
given  by  variational  computations  at  the  normal 
density  of  liquid  helium.   Solid  curve:   Jastrov 
factor   e"^'^'^^^    with   b  =  1.18.   Dotted  line:  hard- 
sphere-like  Jastrow   factor  as  given  by  (9.1) 
and  (9.2)  with   a  =  0. 83. 


Fig.  4     Structure  factor  at  the  normal  density  of 

27  2i 
liquid  helium.   Solid  line:   experiment.    ' 

O0O  :  "reconstructed"  experiment  as  explained 

in  the  text.  *  •  •  :   "exact"  hard  sphere 

structure  factor. 
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This  report  was  prepared  as  an  account  of 
Government  sponsored  work.   Neither  the 
United  States,  nor  the  Commission,  nor  any 
person  acting  on  behalf  of  the  Commission: 

A.  Makes  any  warranty  or  representation, 
express  or  implied,  with  respect  to  the 
accuracy,  completeness,  or  usefulness  of 
the  Information  contained  In  this  report, 
or  that  the  use  of  any  Information, 
apparatus,  method,  or  process  disclosed 
In  this  report  may  not  Infringe  privately 
owned  rights;  or 

B.  Assumes  any  liabilities  with  respect  to 
the  use  of,  or  for  damages  resulting  from 
the  use  of  any  Information,  apparatus, 
method,  or  process  disclosed  In  this 
report. 

As  used  In  the  above,  "person  acting  on  behalf 
of  the  Commission"  Includes  any  employee  or 
contractor  of  the  Commission,  or  employee  of 
such  contractor,  to  the  extent  that  such  em- 
ployee or  contractor  of  the  Commission,  or 
employee  of  such  contractor  prepares,  dis- 
seminates, or  provides  access  to,  any  infor- 
mation pursuant  to  his  employment  or  contract 
with  the  Commission,  or  his  employment  with 
such  contractor. 
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